rm(list = ls())

library(dplyr)
library(magrittr)
library(here)
library(fixest)

# Define relative path to save in the script's directory
file_path <- here()

# You should download "expcor_cls_apsr.rda" to run this script: Tai, Yuehong 'Cassandra'; Hu, Yue; Solt, Frederick, 2022, "Replication Data for: Democracy, Public Support, and Measurement Uncertainty", https://doi.org/10.7910/DVN/XAUF3H, 


load(here('expcor_cls_apsr.rda'))


# Adding V-DEM vars  ----

## Adding vdem variables to further breakdown the liberal component (note that you need "vdemdata" library to run this:)

vdem <- vdemdata::vdem

#e_migdppcln missing
vdemvars <- vdem %>% select(country_name, year, v2x_polyarchy, v2x_liberal, v2xcl_rol, v2x_jucon, v2xlg_legcon, v2clrspct, v2xnp_regcorr,v2x_execorr, v2x_corr, v2xpe_exlecon, v2xpe_exlgender, v2xpe_exlgeo, v2xpe_exlpol, v2xpe_exlsocgr, v2peapsecon, v2peapsgen, v2peapspol, v2peapssoc, v2peasjsoecon, v2pepwrses, v2clacjust, v2x_libdem, v2x_polyarchy, v2jureform, v2jureform_ord, v2jupoatck, v2jureview, e_gdp, e_wbgi_gee, v2exbribe, v2exembez, v2x_pubcorr, e_wbgi_cce, e_ti_cpi, v2x_horacc, v2x_veracc, v2x_diagacc, e_miinflat, v2x_accountability, e_total_fuel_income_pc, e_miurbpop, e_civil_war, e_pt_coup, e_peaveduc, v2x_clphy, v2caviol, e_total_oil_income_pc)

colnames(vdemvars)[which(colnames(vdemvars)=='country_name')] <- 'Country'
colnames(vdemvars)[which(colnames(vdemvars)=='year')] <- 'Year'

## Fixing some country names before merging

vdemvars$Country <- ifelse(vdemvars$Country=='Trinidad and Tobago', 'Trinidad & Tobago', vdemvars$Country)

vdemvars$Country <- ifelse(vdemvars$Country=='United States of America', 'United States', vdemvars$Country)

vdemvars$Country <- ifelse(vdemvars$Country=='Burma/Myanmar', 'Myanmar', vdemvars$Country)

vdemvars$Country <- ifelse(vdemvars$Country=='Palestine/West Bank', 'Palestine', vdemvars$Country)

vdemvars$Country <- ifelse(vdemvars$Country=='Bosnia and Herzegovina', 'Bosnia & Herzegovina', vdemvars$Country)

vdemvars$Country <- ifelse(vdemvars$Country=='Ivory Coast', 'Côte d’Ivoire', vdemvars$Country)

vdemvars$Country <- ifelse(vdemvars$Country=='Türkiye', 'Turkey', vdemvars$Country)


## Add Swiid data (Version 9.7). It is available here: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/LM4OWF

swiid <- read.csv(here('swiid9_7_summary.csv'))

swiid <- swiid %>% select(country, year, gini_disp, gini_disp_se)

# fixing some country names before merging
swiid$country <- ifelse(swiid$country=="Korea", 'South Korea', swiid$country)


swiid$country <- ifelse(swiid$country=='Palestinian Territories', 'Palestine', swiid$country)

swiid$country <- ifelse(swiid$country=='São Tomé and Príncipe', 'Sao Tome and Principe', swiid$country)

swiid$country <- ifelse(swiid$country=='Bosnia and Herzegovina', 'Bosnia & Herzegovina', swiid$country)

swiid$country <- ifelse(swiid$country=="Côte d'Ivoire", 'Côte d’Ivoire', swiid$country)

swiid$country <- ifelse(swiid$country=="Czech Republic", 'Czechia', swiid$country)

swiid$country <- ifelse(swiid$country=="Trinidad and Tobago", 'Trinidad & Tobago', swiid$country)


swiid <- swiid %>% rename(Country = country, Year = year)

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Independent variable: Judicial constraints ----


# Define a function to process each dataset
process_dataset <- function(data, vdemvars, swiid) {
  
  # Merge data
  data <- data %>% 
    rename(Country = country, Year = year) %>%
    left_join(vdemvars, by = c("Country", "Year")) %>%
    left_join(swiid, by = c("Country", "Year")) %>%
    group_by(Country) %>%
    mutate(e_wbgi_cce_imp = zoo::na.approx(e_wbgi_cce, na.rm = FALSE),
        gdp_change_pp = ((e_gdp * 100) / lag(e_gdp)) - 100) %>% 
    mutate(v2x_jucon_lag = lag(v2x_jucon, n = 1),
           gdp_change_pp_lag = lag(gdp_change_pp, n = 1), 
           gini_disp_lag = lag(gini_disp, n = 1),
           e_peaveduc_lag = lag(e_peaveduc, n = 1), 
           e_wbgi_cce_imp_lag = lag(e_wbgi_cce_imp, n = 1),
           v2x_clphy_lag = lag(v2x_clphy, n = 1)) %>% ungroup()
  
  # Define the models
  models <- list(
    feols(SupDem_trim ~ v2x_jucon_lag | Country, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon_lag | Country + Year, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon_lag + gdp_change_pp_lag + gini_disp_lag | Country, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon_lag + gdp_change_pp_lag + gini_disp_lag | Country + Year, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag | Country, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag | Country + Year, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag | Country, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag | Country + Year, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag + v2x_clphy_lag | Country, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag + v2x_clphy_lag | Country + Year, data = data, vcov = ~Country)
  )
  
  # Extract model summaries
  model_summaries <- lapply(models, function(model) {
    # Extract coefficients and model info
    tidy_model <- broom::tidy(model)
    tidy_model$model <- deparse(substitute(model))
    tidy_model
  })
  
  # Combine model summaries into one data frame
  bind_rows(model_summaries, .id = "model_id")
}

# Process all datasets in the list
all_results <- lapply(expcor_cls_apsr, function(dataset) {
  tryCatch(
    process_dataset(dataset, vdemvars, swiid),
    error = function(e) {
      message("Error processing dataset: ", e)
      NULL
    }
  )
})

# Combine all results into a single data frame
final_results <- bind_rows(all_results, .id = "dataset_id")


# Save results for later use
saveRDS(final_results, file = here('regression_results_judicial_const_uncert_lag.rds'))


# Legislative constraints ----

# Define a function to process each dataset
process_dataset_leg <- function(data, vdemvars, swiid) {
  
  # Merge data
  data <- data %>% 
    rename(Country = country, Year = year) %>%
    left_join(vdemvars, by = c("Country", "Year")) %>%
    left_join(swiid, by = c("Country", "Year")) %>%
    group_by(Country) %>%
    mutate(
      e_wbgi_cce_imp = zoo::na.approx(e_wbgi_cce, na.rm = FALSE),
      gdp_change_pp = ((e_gdp * 100) / lag(e_gdp)) - 100
    ) %>% mutate(v2xlg_legcon_lag = lag(v2xlg_legcon, n = 1),
                 gdp_change_pp_lag = lag(gdp_change_pp, n = 1), 
                 gini_disp_lag = lag(gini_disp, n = 1),
                 e_peaveduc_lag = lag(e_peaveduc, n = 1), 
                 e_wbgi_cce_imp_lag = lag(e_wbgi_cce_imp, n = 1),
                 v2x_clphy_lag = lag(v2x_clphy, n = 1)) %>% 
    ungroup()
  
  # Define the models
  models <- list(
    feols(SupDem_trim ~ v2xlg_legcon_lag | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon_lag | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon_lag + gdp_change_pp_lag + gini_disp_lag | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon_lag + gdp_change_pp_lag + gini_disp_lag | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag + v2x_clphy_lag | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag + v2x_clphy_lag | Country + Year, data = data, vcov = ~Country)
    
  )
  
  # Extract model summaries
  model_summaries <- lapply(models, function(model) {
    # Extract coefficients and model info
    tidy_model <- broom::tidy(model)
    tidy_model$model <- deparse(substitute(model))
    tidy_model
  })
  
  # Combine model summaries into one data frame
  bind_rows(model_summaries, .id = "model_id")
}

# Process all datasets in the list
all_results_leg <- lapply(expcor_cls_apsr, function(dataset) {
  tryCatch(
    process_dataset_leg(dataset, vdemvars, swiid),
    error = function(e) {
      message("Error processing dataset: ", e)
      NULL
    }
  )
})

# Combine all results into a single data frame
final_results_leg <- bind_rows(all_results_leg, .id = "dataset_id")

# Save results for later use
saveRDS(final_results_leg, file = here('regression_results_legislative_const_uncert_lag.rds'))

# Ind liberties ----

# Define a function to process each dataset
process_dataset_indlib <- function(data, vdemvars, swiid) {
  
  # Merge data
  data <- data %>% 
    rename(Country = country, Year = year) %>%
    left_join(vdemvars, by = c("Country", "Year")) %>%
    left_join(swiid, by = c("Country", "Year")) %>%
group_by(Country) %>%
    mutate(
      e_wbgi_cce_imp = zoo::na.approx(e_wbgi_cce, na.rm = FALSE),
      gdp_change_pp = ((e_gdp * 100) / lag(e_gdp)) - 100
    ) %>% 
    mutate(v2xcl_rol_lag = lag(v2xcl_rol, n = 1),
           gdp_change_pp_lag = lag(gdp_change_pp, n = 1), 
           gini_disp_lag = lag(gini_disp, n = 1),
           e_peaveduc_lag = lag(e_peaveduc, n = 1), 
           e_wbgi_cce_imp_lag = lag(e_wbgi_cce_imp, n = 1),
           v2x_clphy_lag = lag(v2x_clphy, n = 1)) %>% ungroup()
  
  # Define the models
  models <- list(
    feols(SupDem_trim ~ v2xcl_rol_lag | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol_lag | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol_lag + gdp_change_pp_lag + gini_disp_lag | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol_lag + gdp_change_pp_lag + gini_disp_lag | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag + v2x_clphy_lag | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag + v2x_clphy_lag | Country + Year, data = data, vcov = ~Country)
    
  )
  
  # Extract model summaries
  model_summaries <- lapply(models, function(model) {
    # Extract coefficients and model info
    tidy_model <- broom::tidy(model)
    tidy_model$model <- deparse(substitute(model))
    tidy_model
  })
  
  # Combine model summaries into one data frame
  bind_rows(model_summaries, .id = "model_id")
}

# Process all datasets in the list
all_results_indlib <- lapply(expcor_cls_apsr, function(dataset) {
  tryCatch(
    process_dataset_indlib(dataset, vdemvars, swiid),
    error = function(e) {
      message("Error processing dataset: ", e)
      NULL
    }
  )
})

# Combine all results into a single data frame
final_results_indlib <- bind_rows(all_results_indlib, .id = "dataset_id")

# Save results for later use
saveRDS(final_results_indlib, file = here('regression_results_indliberties_const_uncert_lag.rds'))


